Granger Causality via Vector Auto-Regression Tuned for FMRI Data Analysis

Introduction: In addition to localizing task-related cortical and subcortical activity, FMRI data can be used to model various forms of interactions among active regions of interest (ROIs) in the brain. Granger causality (GC) is a connectivity modeling approach with the following premises: (a) a region causally connected to another would leave in the latter a time-lagged version of its own signal components; and (b) the latency in the signal indicates a causal relation. In a typical multivariate GC analysis, the causality among the pre-selected ROIs is determined in a vector auto-regressive (VAR) model with appropriate lags (time delay order) over which connectivity is to be examined. A second-stage analysis can then be conducted to examine group differences in connectivity. We present a multivariate implementation of GC that takes into account the idiosyncrasies of FMRI data and provides for inference on the sign of correlation, evolution of the network over multiple lags, model order tuning, and model diagnosis. We also propose a valid group analysis per lag based on signed path coefficients to reveal the network at the group level. Methods: We start with an extended VAR(p) model Y(t)=α+A1Y(t-1)+...+ApY(t-p)+BZ(t)+ε(t), where: Y(t) is an n×1 vector containing the BOLD response of n ROIs at time t; Z(t) is an q×1 vector containing exogenous variables (covariates or confounds) at time t; p and q are the number of lags and confounds respectively; α is an n×1 vector modeling the baseline; Ai is an n×n matrix with each element reflecting the effect from one ROI to another at time lag i (i=1,...,p); B is an n×q matrix with each element modeling the effect of a confound on an ROI; and ε(t) is the error term. The confounds can be baseline drift, head motion parameters, tasks of no interest, physiological measurements, or dummy variables corresponding to time breaks in the time series due to multiple blocks/runs/sessions and to time points censored because of signal irregularities. Unlike previous implementations [1, 2] of the standard VAR(p) model Y(t)=α+A1Y(t-1)+...+ApY(t-p)+ε(t) with significance based on Fstatistics across multiple lags, our approach with an extended model allows us to perform whole network causality analysis while simultaneously modeling confounding effects, instead of regressing them out prior to the GC analysis. Furthermore, our implementation carries correlation signs at different lags all the way to the group level analysis, because reducing correlations to F-statistics and collapsing networks across lags can significantly alter the identified network. Four criteria (AIC, FPE, HQ, and SC) are provided to assist the user in selecting appropriate lag levels (order) of the model, and one network per lag is identified based on path coefficients and their t-statistics, potentially showing the network evolution through time (lags). Lastly, a set of model diagnosis tests are available to check for stationarity (testing roots of characteristic polynomial), residual normality tests (Gaussian process, skewness, and kurtosis), residual autocorrelation (portmanteau, Breusch-Godfrey, and Edgerton-Shukur tests), autoregressive conditional heteroskedasticity (ARCH) test for time-varying volatility, and structure stability (or stationarity detection). Our modeling approach was validated at the individual subject level by comparing with [2] when using the same model and input. To illustrate the importance of preserving the correlation sign, we modeled the network among three ROIs—dorsolateral prefrontal cortex (DLPFC), dorsolateral caudate, and thalamus—that comprise the cognitive division of the cortico-striatal-thalamic loop, a well-validated neural circuit [3]. We collected 250 volume images during eye-closed rest [cf. ref. 4: 18 axial slices, spiral in-out; TR=1.2s, TE=30ms, FOV=22cm, matrix=64×64, slice thickness=5mm, flip angle=77°], from each of 20 participants. Results: For this analysis, we adopted a lag of 1, as suggested by the model tuning step of our software written in R [5, 6]. The left half of Figure 1 shows the group network obtained with a one-sample t-test using log-transformed F-statistics of individual subjects for each connection, as typically practiced in the literature [e.g., 1]. The network in the right half of Figure 1 was obtained with a one-sample t-test using the signed path coefficients instead. Note the marked differences between the two networks: none of the inter-regional connections were identified with transformed F-statistics as input, and the information about the correlation sign was totally lost, which can be attributed to either ignoring the path coefficient sign, the arbitrary transformation of the F-statistics, or both.